---
title: "Replication Materials for Unveiled"
author: "Aengus Bridgman, Costin Ciobanu, Aaron Erlich"
date: "30/01/2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

# Data and packages

```{r}

library(zoo)
library(MASS)
library(texreg)
library(miceadds)
library(matchMulti)
library(robustbase)
library(estimatr)
library(gridExtra)
library(cobalt)
library(konfound)
library(rdrobust)
library(tidyverse)
library(tidyquant)
library(fastDummies)

options(digits = 2, scipen=10)

theme <- theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.position="top",
    strip.background = element_blank(),
    panel.background=element_blank(),
    panel.grid.minor=element_blank(),
    panel.grid.major= element_line(color = "seashell1"),
    plot.background=element_blank()
  )

grobtheme <- ttheme_minimal(
    core = list(fg_params=list(cex = 0.5)),
    colhead = list(fg_params=list(cex = 0.5)),
    rowhead = list(fg_params=list(cex = 0.5))
    )


dir.create("outputs/")
dir.create("outputs/tables")
dir.create("outputs/figures")

```


```{r}

CES <- readRDS(file = "data/ready/CES.rds") %>% filter(discard != 1)
LPP  <- readRDS(file = "data/ready/LPP.rds")
media <- readRDS(file = "data/ready/media.rds")
media_sentiment = readRDS("data/ready/media_sentiment.rds")

```

# MAIN BODY

## LPP Quebec Stat
```{r}

LPP %>%
  filter(vote_for_all_parties != "undecided") %>% 
  group_by(court, quebec) %>% 
  summarize(ndp_vote = mean(ndp_vote, na.rm = TRUE), n = n()) 

```

## fig1
```{r}

fig1 <- media %>%
  dplyr::select(date, quebec, economy, niqab) %>%
  rename(`Niqab coverage` = niqab,
         `Economy coverage` = economy) %>%
  ungroup() %>%
  gather(key = "key", value = "value", -quebec, -date) %>%
  mutate(quebec = ifelse(quebec == 1, "Quebec","Rest of Canada"))

p = ggplot(fig1, aes(x = date, y = value, linetype = key, col = key)) +
  scale_x_date(breaks = seq(as.Date("2015-07-05"),as.Date("2015-10-19"),7),
               limits = c(as.Date("2015-07-05"),as.Date("2015-10-19")),
               date_labels = "%b %e") +
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-09-15"))),
             linetype = 2, size = 0.8, col = "darkgrey") +
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-08-02"))),
             linetype = 2, size = 0.8, col = "darkgrey") +
  facet_grid(quebec~., scales = "free_y") +
  geom_ma(aes(lty = key, col = key), ma_fun = SMA, n = 7, cex = 1) +
  scale_linetype_manual(name = "", values = c("solid","solid")) +
  theme + 
  labs(y = "Rolling number of articles per newspaper \n in Quebec and the Rest of Canada", x = "") +
  scale_color_manual(name = "", values = c("darkgrey","black")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
   theme(legend.key.size = grid::unit(5, "lines"))

dat_text <- data.frame(key = c("Economy coverage","Economy coverage"),
                       value = c(4,4),
                       quebec = c("Quebec","Quebec"),
                       label = c("Campaign begins","Court Ruling"),
                       xpos = c(as.Date("2015-07-25"),as.Date("2015-09-09")))

p + geom_text(data = dat_text,
              mapping = aes(x = xpos, y = value, label = label),
              size = 2,
              col = "black",
            show.legend=FALSE) +
  coord_cartesian(ylim=c(0,4))

ggsave("outputs/figures/fig1.png", width = 6.5, height = 5)

```

## tab1
```{r}

tab1 <- list()

tab1[[1]] <- lm_robust(data = filter(CES, web == 1, !is.na(region)) %>% mutate(ndp_vote = ndp_vote*100),
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + french + age + region_AT + region_PR + region_BC + 
                         working + student + retired + education +
                         vote_2011_ndp + court*quebec)

tab1[[2]] <- lm_robust(data = filter(CES, web == 1, !is.na(region)) %>% mutate(ndp_vote = ndp_vote*100), 
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + age + french + region_AT + region_PR + region_BC + 
                         working + student + retired + education + 
                         vote_2011_ndp + court*quebec + court_linear*quebec)

tab1[[3]] <-  lm.cluster(data = filter(CES, web == 1, !is.na(region)) %>% mutate(ndp_vote = ndp_vote*100), 
                         cluster = "date",
                         formula = ndp_vote ~ female + age + french + region_AT + region_PR + region_BC +
                           working + student + retired + education +
                           vote_2011_ndp + niqab_7*quebec)

screenreg(list(tab1[[1]],tab1[[2]],tab1[[3]]$lm_res), include.ci = FALSE,
          omit.coef = "(female)|(age)|(french)|(region)|
          (PR)|(working)|(student)|(retired)|(education)")

# ERROR: 1 coefficient  not defined because the design matrix is rank deficient
# This is because both region and an interaction with Quebec are included in the model
texreg(list(tab1[[1]],tab1[[2]],tab1[[3]]$lm_res),
       scalebox = 1, float.pos = "ht!", label = "tab1", include.ci = FALSE,
       dcolumn = TRUE,  use.packages = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP", 
       caption.above = TRUE, single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)",
       custom.note = ("\\parbox{0.9\\linewidth}{\\vspace{2pt}%stars. Linear probability models for DID estimations with robust standard errors for Models 1 and 2 and clustered standard errors for Model 3 in parentheses. Dependent variable is vote intention for the NDP (binary variable). All models use full web sample.}"),
       custom.model.names = c("1: Binary DID","2: Linear Trend ","3: 7-day media"),
       custom.coef.names = c("Constant","Voted NDP 2011",
                             "Ruling","Quebec","Ruling x Quebec","Trend","Trend x Quebec",
                             "7-day niqab","7-day niqab x Quebec"),
       reorder.coef = c(5,7,9,1,2,3,4,6,8),
       groups = list("DID effects" = 1:3, "Other coefficients" = 4:9),
       file = "outputs/tables/tab1.tex")


```

## Media sentiment
```{r}

daily_sentiment = media_sentiment %>% 
  group_by(date) %>% 
  summarize(negative_lpc = sum(negative_lpc, na.rm = TRUE),
            positive_lpc = sum(positive_lpc, na.rm = TRUE),
            negative_ndp = sum(negative_ndp, na.rm = TRUE),
            positive_ndp = sum(positive_ndp, na.rm = TRUE)) %>% 
  ungroup() %>% 
  mutate(prop_pos_lpc = positive_lpc/(positive_lpc+negative_lpc),
         prop_neg_lpc = negative_lpc/(positive_lpc+negative_lpc),
         prop_pos_ndp = positive_ndp/(positive_ndp+negative_ndp),
         prop_neg_ndp = negative_ndp/(positive_ndp+negative_ndp),
         lib_net_tone = log((prop_pos_lpc+0.05)/(prop_neg_lpc+0.05)),
         ndp_net_tone = log((prop_pos_ndp+0.05)/(prop_neg_ndp+0.05)))

mean(daily_sentiment$prop_neg_lpc, na.rm = T)
mean(daily_sentiment$prop_neg_ndp, na.rm = T)
mean(media_sentiment$prop_negative_lpc, na.rm = T)
mean(media_sentiment$prop_negative_ndp, na.rm = T)

# T-test to determine the statistical difference over the period
t.test(media_sentiment$prop_positive_lpc, media_sentiment$prop_positive_ndp) # each article
t.test(daily_sentiment$lib_net_tone, daily_sentiment$ndp_net_tone) # day aggregate

# T-test to determine the statistical difference over the period (0.004)
t.test(media_sentiment$lib_net_tone, media_sentiment$ndp_net_tone) # each article
t.test(daily_sentiment$lib_net_tone, daily_sentiment$ndp_net_tone) # day aggregate

```

## Media mentions
```{r}

# T-test to determine the statistical difference over the period 
t.test(media_sentiment$lpc_mentions, media_sentiment$ndp_mentions)

```

## Panel
```{r}

filter(CES, !is.na(p_ban_niqab), quebec == 1) %>%
  dplyr::select(quebec, court, p_ban_niqab,ndp_vote) %>%
  group_by(quebec, court, p_ban_niqab) %>%
  summarize(mean = mean(ndp_vote, na.rm = TRUE), n = n())

t.test(filter(CES, court == 0, p_ban_niqab == 0, quebec == 1)$ndp_vote,
       filter(CES, court == 0, p_ban_niqab == 1, quebec == 1)$ndp_vote)

t.test(filter(CES, court == 1, p_ban_niqab == 0, quebec == 1)$ndp_vote,
       filter(CES, court == 1, p_ban_niqab == 1, quebec == 1)$ndp_vote)

t.test(filter(CES, court == 0, quebec == 1)$ndp_voted,
       filter(CES, court == 1, quebec == 1)$ndp_voted)

```

# APPENDIX

## A-1 Data sets employed - information and descriptive statistics

### tabA1
```{r}

CES_descriptives <- filter(CES, web == 1) %>% 
  dplyr::select(ndp_vote, lib_vote, con_vote,
                p_ban_niqab,
                age, female, french, quebec,
                working, student, retired, 
                no_hs, hs, bach, grad, 
                vote_2011_ndp,
                feel_con, feel_lib, feel_ndp, feel_blq, 
                feel_hrpr, feel_trud, feel_mulc, feel_ducp,
                court, court_linear) %>%
  psych::describe(skew = TRUE, check = TRUE) %>%
  dplyr::select(n, mean, sd, median, min, max) %>%
  as_tibble() %>%
  rownames_to_column(.)

colnames(CES_descriptives) <- c("Variable","n","Mean","SD","Median","Min","Max")

CES_descriptives[,'Variable'] <- c(
  "NDP Vote Intention",
  "Liberal Vote Intention",
  "Conservative Vote Intention",
  "Favor of Niqab Ban",
  "Age","Female","French", "Quebec",
  "Working","Student","Retired",
  "No High School","High School",
  "Bachelor's Degree","Graduate Studies",
  "Vote for NDP in 2011",
  "Feeling towards Conservative Party",
  "Feeling towards Liberal Party",
  "Feeling towards NDP",
  "Feeling towards the Bloc Quebecois",
  "Feeling towards Harper",
  "Feeling towards Trudeau",
  "Feeling towards Mulcair",
  "Feeling towards Duceppe",
  "Post-Court decision",
  "Post-Court decision linear trend"
)

print(xtable::xtable(CES_descriptives, 
                     digits = c(0,0,0,2,2,2,0,0),
                     label = "tabA1",
                     caption = "Descriptive statistics for CES data (web sample)"),
      file="outputs/tables/tabA1.tex",
      caption.placement = "top", scalebox='0.9')

```

### tabA2
```{r}

media_descriptives <- filter(media, quebec == 0) %>%
  ungroup() %>%
  dplyr::select(niqab_num, niqab_7, cum_niqab, econo_num, econo_num_7, cum_econo) %>%
  rename(`Daily niqab coverage` = niqab_num,
         `Rolling sum of niqab coverage (7 days)` = niqab_7,
         `Cumulative coverage of niqab` = cum_niqab,
         `Daily economy coverage` = econo_num,
         `Rolling sum of economy coverage (7 days)` = econo_num_7,
         `Cumulative coverage of economy` = cum_econo) %>%
  psych::describe(skew = TRUE, check = TRUE) %>%
  dplyr::select(n, mean, sd, median, min, max) %>%
  as_tibble() %>%
  rownames_to_column(.)

colnames(media_descriptives) <- c("Variable","n","Mean","SD","Median","Min","Max")

print(xtable::xtable(media_descriptives,
                     digits = c(0,0,0,2,2,2,0,0),
                     label = "tabA2",
                     caption = "Descriptive statistics for media data (rest of Canada)"),
      file="outputs/tables/tabA2.tex", 
      caption.placement = "top", scalebox='0.9')

```

### tabA3
```{r}

media_descriptives <- filter(media, quebec == 1) %>%
  ungroup() %>%
  dplyr::select(niqab_num, niqab_7, cum_niqab, econo_num, econo_num_7, cum_econo) %>%
  rename(`Daily niqab coverage` = niqab_num,
         `Rolling sum of niqab coverage (7 days)` = niqab_7,
         `Cumulative coverage of niqab` = cum_niqab,
         `Daily economy coverage` = econo_num,
         `Rolling sum of economy coverage (7 days)` = econo_num_7,
         `Cumulative coverage of economy` = cum_econo) %>%
  psych::describe(skew = TRUE, check = TRUE) %>%
  dplyr::select(n, mean, sd, median, min, max) %>%
  as_tibble() %>%
  rownames_to_column(.)

colnames(media_descriptives) <- c("Variable","n","Mean","SD","Median","Min","Max")

print(xtable::xtable(media_descriptives,
                     digits = c(0,0,0,2,2,2,0,0),
                     label = "tabA3",
                     caption = "Descriptive statistics for media data (Quebec-only)"), file="outputs/tables/tabA3.tex",
      caption.placement = "top", scalebox='0.9')

```

## A-2 Evolution of vote intention and party identification

### figA1
```{r}


ggdat <- filter(CES, date <= "2015-10-18" & !is.na(date)) %>%
  filter(!is.na(quebec)) %>%
  group_by(date) %>%
  summarise(`CPC Vote` = mean(con_vote, na.rm = TRUE),
            `NDP Vote` = mean(ndp_vote, na.rm = TRUE),
            `LPC Vote` = mean(lib_vote, na.rm = TRUE),
            `CPC PID` = mean(con_pid, na.rm = TRUE),
            `NDP PID` = mean(ndp_pid, na.rm = TRUE),
            `LPC PID` = mean(lib_pid, na.rm = TRUE)) %>%
  gather(key, value, -date) %>%
  mutate(pid = ifelse(key %in% c("CPC PID","LPC PID","NDP PID"),"PID","Vote intention"),
         party = case_when(
           key %in% c("CPC PID","CPC Vote") ~ "CPC",
           key %in% c("LPC PID","LPC Vote") ~ "LPC",
           key %in% c("NDP PID","NDP Vote") ~ "NDP"
         ))

ggdat %>%
  ggplot(., aes(x = date, y = value, alpha = pid)) +
  geom_ma(aes(col = key), ma_fun = SMA, n = 7, lty = 1, cex = 2) +
  theme + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "", y = "% Respondents") + 
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-09-15"))), linetype = 3) +
  scale_color_manual(name = "", guide = FALSE, values = c("#1A4782","#1A4782",
                                           "#D71920","#D71920",
                                           "#F37021","#F37021")) +
  scale_alpha_manual(name = "", values=c(0.5, 1)) +
  facet_grid(.~party)

ggsave("outputs/figures/figA1.png", width = 7.5, height = 4)


```

### figA2
```{r}

ggdat <- filter(CES, date <= "2015-10-18" & !is.na(date)) %>%
  filter(!is.na(quebec)) %>%
  group_by(date, quebec) %>%
  summarise(`CPC Vote` = mean(con_vote, na.rm = TRUE),
            `NDP Vote` = mean(ndp_vote, na.rm = TRUE),
            `LPC Vote` = mean(lib_vote, na.rm = TRUE),
            `BQ Vote` = mean(blq_vote, na.rm = TRUE),
            `CPC PID` = mean(con_pid, na.rm = TRUE),
            `NDP PID` = mean(ndp_pid, na.rm = TRUE),
            `LPC PID` = mean(lib_pid, na.rm = TRUE),
            `BQ PID` = mean(blq_pid, na.rm = TRUE)) %>%
  gather(key, value, -date, -quebec) %>%
  mutate(quebec = ifelse(quebec == 1, "Quebec", "ROC"),
         pid = ifelse(key %in% c("CPC PID","LPC PID","NDP PID","BQ PID"),"PID","Vote intention"),
         party = case_when(
           key %in% c("CPC PID","CPC Vote") ~ "CPC",
           key %in% c("LPC PID","LPC Vote") ~ "LPC",
           key %in% c("NDP PID","NDP Vote") ~ "NDP",
           key %in% c("BQ PID","BQ Vote") ~ "BQ"
         ))

filter(ggdat, quebec == "Quebec") %>%
  ggplot(., aes(x = date, y = value, alpha = pid)) +
  geom_ma(aes(col = key), ma_fun = SMA, n = 7, lty = 1, cex = 2) +
  theme + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "", y = "% Respondents") + 
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-09-15"))), linetype = 3) +
  scale_color_manual(name = "", guide = FALSE, values = c("#33B2CC","#33B2CC",
                                                          "#1A4782","#1A4782",
                                                          "#D71920","#D71920",
                                                          "#F37021","#F37021")) +
  scale_alpha_manual(name = "", values=c(0.5, 1)) +
  facet_grid(.~party)

ggsave("outputs/figures/figA2.png", width = 7.5, height = 4)

```

### figA3
```{r}

ggdat <- filter(CES, date <= "2015-10-18" & !is.na(date), quebec == 1) %>%
  filter(!is.na(quebec)) %>%
  group_by(date) %>%
  summarise(`CPC Vote` = mean(vote_2011_con, na.rm = TRUE),
            `NDP Vote` = mean(vote_2011_ndp, na.rm = TRUE),
            `LPC Vote` = mean(vote_2011_lib, na.rm = TRUE)) %>%
  gather(key, value, -date) %>%
  mutate(pid = ifelse(key %in% c("CPC PID","LPC PID","NDP PID"),"PID","Vote intention"),
         party = case_when(
           key %in% c("CPC PID","CPC Vote") ~ "CPC",
           key %in% c("LPC PID","LPC Vote") ~ "LPC",
           key %in% c("NDP PID","NDP Vote") ~ "NDP"
         ))

filter(ggdat) %>%
  ggplot(., aes(x = date, y = value)) +
  geom_ma(aes(col = key), ma_fun = SMA, n = 7, lty = 1, cex = 2) +
  theme + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "", y = "% Respondents") + 
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-09-15"))), linetype = 3) +
  scale_color_manual(name = "", guide = FALSE, values = c("#1A4782","#D71920","#F37021")) +
  facet_grid(.~party)

ggsave("outputs/figures/figA3.png", width = 7.5, height = 4)

```

## B-1 The context of the 2015 Canadian federal election

### figB1
```{r}

ggdat <- filter(CES, date <= "2015-10-18" & !is.na(date)) %>%
  dplyr::select(con_vote, lib_vote, ndp_vote, date) %>%
  group_by(date) %>%
  summarise(`CPC Vote` = mean(con_vote, na.rm = TRUE),
            `NDP Vote` = mean(ndp_vote, na.rm = TRUE),
            `LPC Vote` = mean(lib_vote, na.rm = TRUE)) %>%
  gather(key, value, -date) %>%
  mutate(value = value*100)
  
ggplot(ggdat, aes(x = date, y = value)) +
  geom_ma(aes(col = key), ma_fun = SMA, n = 7, lty = 1, cex = 2) +
  theme + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "", y = "% Respondents") + 
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-09-15"))), linetype = 4) +
  annotate("text", x = as.Date("2015-09-22"), y = 33, label = "Court Decision") + 
  scale_color_manual(name = "",values = c("#1A4782","#D71920","#F37021"))

ggsave("outputs/figures/figB1.png", width = 7.5, height = 4)

```

### figB2
```{r}

ggdat <- filter(CES, date <= "2015-10-18" & !is.na(date)) %>%
  dplyr::select(con_vote, lib_vote, ndp_vote, blq_vote, date, quebec) %>%
  filter(!is.na(quebec)) %>%
  group_by(date, quebec) %>%
  summarise(`CPC Vote` = mean(con_vote, na.rm = TRUE),
            `NDP Vote` = mean(ndp_vote, na.rm = TRUE),
            `LPC Vote` = mean(lib_vote, na.rm = TRUE),
            `BLQ Vote` = mean(blq_vote, na.rm = TRUE)) %>%
  gather(key, value, -date, -quebec) %>%
  mutate(value = value*100) %>%
  mutate(qc = ifelse(quebec == 1, "Quebec", "ROC"))
  
filter(ggdat, qc == "Quebec") %>%
  ggplot(aes(x = date, y = value)) +
  geom_ma(aes(col = key), ma_fun = SMA, n = 7, lty = 1, cex = 2) +
  theme + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "", y = "% Respondents") + 
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-09-15"))), linetype = 3) +
  annotate("text", x = as.Date("2015-09-22"), y = 38, label = "Court Decision") + 
  scale_color_manual(name = "",values = c("#33B2CC","#1A4782","#D71920","#F37021"))

ggsave("outputs/figures/figB2.png", width = 7.5, height = 4)

```

## B-2 Religious symbols and the Quebec distinctiveness

### figB3
```{r}

CES %>%
  group_by(quebec) %>%
  dplyr::select(quebec,p_ban_niqab,p_musl) %>%
  mutate(p_ban_niqab = p_ban_niqab, p_musl = p_musl/100) %>%
  rename(`Favour banning niqab` = p_ban_niqab,
         `Feel about Muslims` = p_musl) %>%
  summarise_all(funs(mean = mean(., na.rm = TRUE), sd = sd(., na.rm = TRUE), n = sum(!is.na(.)))) %>%
  gather(key= "key", value = "value", -quebec) %>%
  separate(key, into = c("variable","operation"), sep = "_") %>%
  spread(key = "operation", value = "value") %>%
  mutate(se=sd/sqrt(n))  %>%
  mutate(ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  mutate(quebec = ifelse(quebec == 1, "Quebec","Rest of Canada")) %>%
  ggplot(aes(fill = quebec)) +
  geom_bar(aes(x=quebec, y=mean), stat="identity", alpha=0.5, width = 0.5, color = "black") +
  geom_errorbar(aes(x=quebec, ymin=mean-ic, ymax=mean+ic),
                colour="black", alpha=0.9, size=1, width = 0) +
  facet_wrap(variable~.) +
  theme +
  scale_fill_manual(values = c("slategrey","lightgrey")) +
  labs(x = "", y = "Variables rescaled for comparability (0-1)") +
  guides(fill = FALSE)

ggsave("outputs/figures/figB3.png", width = 6.5, height = 5)

```

## B-3 Alternative explanations

### Descriptives
```{r}

### senate scandal
CES %>%
  group_by(quebec) %>%
    summarize(hear_senate = mean(isslist_sena, na.rm = TRUE),
              care_senate = mean(isscare_sena, na.rm = TRUE),
              n = n())

### budget
CES %>%
  mutate(deficit = ifelse(deficit == 2, 1,0)) %>%
  group_by(quebec) %>%
  summarize(deficit = mean(deficit, na.rm = TRUE))

t.test(filter(CES %>% mutate(deficit = ifelse(deficit == 2, 1,0)), quebec == 0, web == 1)$deficit,
       filter(CES %>% mutate(deficit = ifelse(deficit == 2, 1,0)), quebec == 1, web == 1)$deficit)

budget <- list()

budget[[1]] <- lm(data = filter(CES %>% mutate(deficit = ifelse(deficit == 2, 1,0)), web == 1),
                100*ndp_vote ~ female + age + french + deficit + 
                  region_AT + region_PR + region_BC +working + student + retired + education +
                  vote_2011_ndp + court*quebec)
budget[[2]] <- lm(data = filter(CES %>% mutate(deficit = ifelse(deficit == 2, 1,0)), web == 1),
                100*ndp_vote ~ female + age + french + deficit + 
                  region_AT + region_PR + region_BC +working + student + retired + education + 
                  vote_2011_ndp + court*quebec + court_linear*quebec)

screenreg(budget)

### leader evaluations
evals <- CES %>% 
  filter(date >= "2015-10-04", date <= "2015-10-18", !is.na(date), quebec == 1) %>%
  group_by(date) %>%
  summarise(feel_mulc = mean(feel_mulc, na.rm = TRUE),
            feel_trud = mean(feel_trud, na.rm = TRUE))

t.test(evals$feel_mulc, evals$feel_trud)

evals <- CES %>% 
  filter(date >= "2015-10-04", date <= "2015-10-18", !is.na(date), quebec == 1)

t.test(evals$feel_mulc, evals$feel_trud)

### immigration
t.test(filter(CES, quebec == 0, web == 1)$isscare_immg,
       filter(CES, quebec == 1, web == 1)$isscare_immg)

```

### tabB1
```{r}

CES_strategic = CES %>% 
  filter(web == 1, date <= "2015-10-12") %>% 
  mutate(ndp_vote = 100*ndp_vote)

strategic <- list()

strategic[[1]] <- lm_robust(data = CES_strategic, 
                            ndp_vote ~ female + age + french + 
                              region_AT + region_PR + region_BC +working + student + retired + education +
                              vote_2011_ndp + court*quebec)
strategic[[2]] <- lm_robust(data = CES_strategic, 
                     ndp_vote ~ female + age + french + 
                       region_AT + region_PR + region_BC +working + student + retired + education + 
                       vote_2011_ndp + court_linear*quebec)
strategic[[3]] <- lm_robust(data = CES_strategic, 
                     ndp_vote ~ female + age + french + 
                       region_AT + region_PR + region_BC +working + student + retired + education + 
                       vote_2011_ndp + court*quebec + court_linear*quebec)

screenreg(list(strategic[[1]],strategic[[2]],strategic[[3]]),
          omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)", 
          include.ci = FALSE)

texreg(strategic, scalebox = 1, float.pos = "ht!", label = "tabB1", 
       dcolumn = TRUE,  use.packages = FALSE, include.ci = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP", 
       caption.above = TRUE, single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)",
       custom.note = ("\\parbox{0.9\\linewidth}{\\vspace{2pt}%stars. Linear probability models for DID estimations with robust standard errors in parentheses. Dependent variable is vote intention for the NDP in the 2015 Canadian federal election (binary). All models use web sample but exclude the last week of the campaign.}"),
       custom.model.names = c("1: DID Quebec","2: Trend Quebec","3: Both"),
       custom.coef.names = c("Constant","Voted NDP 2011", "Ruling","Quebec",
                             "Ruling x Quebec","Trend","Trend x Quebec","Trend x Quebec"),
       reorder.coef = c(1,2,3,4,6,5,7),
       groups = list(" " = 1:5,"DID coefficients" = 6:7),
       file = "outputs/tables/tabB1.tex")

```

## C-1 Exogeneity of the court ruling

### figC1
```{r}

google_trend <- read.csv(file = "data/raw/google_trend.csv")

google_trend <- google_trend %>% 
  mutate(date = as.Date(Day),
         proportion = niqab)

ggplot(google_trend, aes(x = date, y = proportion)) +
  geom_ma(col = "black", size = 1, n = 7, linetype = "solid") +
  theme +
  scale_x_date(breaks = seq(as.Date("2015-07-05"),as.Date("2015-10-19"),7),
               limits = c(as.Date("2015-07-05"),as.Date("2015-10-19")),
               date_labels = "%b %e") +
  geom_vline(aes(xintercept = as.numeric(as.Date("2015-09-15"))),
             linetype = 2, size = 0.8, col = "darkgrey") +
  labs(y = "Rolling proportion of maximum number of \n Google searches over time period under \n examination (n = 7)", x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    annotate("text", x = as.Date("2015-09-01"), y = 60, label = "Court Decision")

ggsave("outputs/figures/figC1.png", width = 7, height = 5)

```

## C-2 Causal effect of court ruling on niqab

### tabC1
```{r}

daily_media <- media %>%
  dplyr::select(date, niqab_p, econo_p, quebec) %>%
  gather(key = "key", value = "value", -date, -quebec) %>%
  mutate(court = ifelse(date > "2015-09-15", 1, 0),
         court_linear = ifelse(date > "2015-09-15", as.numeric(date - 16692), 0),
         day_of_campaign = (as.numeric(date) - as.numeric(as.Date("2015-08-01"),
                                                          na.rm = T)) + 1) %>%
  filter(date > as.Date("2015-08-01"), date < as.Date("2015-10-19"))

media_models <- list()

media_models[[1]] <- lm_robust(data = filter(daily_media, quebec == 1) %>% mutate(value = value*100), 
                               ci = TRUE, se_type = "HC2",
                               value ~ day_of_campaign + key*court + key*court_linear)
media_models[[2]] <- lm_robust(data = filter(daily_media, quebec == 0) %>% mutate(value = value*100),
                               ci = TRUE, se_type = "HC2",
                               value ~ day_of_campaign + key*court + key*court_linear)

screenreg(media_models, single.row = TRUE, include.ci = FALSE)

texreg(media_models, scalebox = 1, float.pos = "ht!", label = "tab:media_did", 
       dcolumn = TRUE,  use.packages = FALSE, include.ci = FALSE,
       caption = "Niqab and economy coverage during election", caption.above = TRUE,
       custom.note = ("\\parbox{0.75\\linewidth}{\\vspace{2pt}%stars. OLS models for DID estimations with robust standard errors in parentheses. Dependent variable is number of articles published per day (continuous variable).}"),
       single.row = TRUE,
       custom.model.names = c("Quebec", "Rest of Canada"),
       custom.coef.names = c("Constant", "Daily campaign trend","Niqab Coverage",
                             "Ruling","Trend",
                             "Niqab x Ruling","Niqab x Trend"),
       groups = list(" " = 1:5,"DID coefficients" = 6:7),
       file = "outputs/tables/tabC1.tex")

```

## D-2 Balance tests

### figD1
```{r}

temp <- filter(CES, web == 1) %>%
  dplyr::select(web, court, age, female, french, region, 
         working, retired, student, 
         vote_2011_ndp, #ndp_voted, ndp_pid, lib_pid, lib_vote,
         education, income_full, religion2, p_ban_niqab,
         interest, knowledge) #ad_exposure

temp %>% 
  bal.tab(., treat = .$court) %>%
  .[[1]] %>%
  rownames_to_column() %>%
  filter(rowname != "court") %>%
  mutate(rowname = case_when(
    rowname == "age" ~ "Age",
    rowname == "working" ~ "Employment (Working)",
    rowname == "vote_2011_ndp" ~ "Voted for NDP in 2011",
    rowname == "student" ~ "Employment (Student)",
    rowname == "retired" ~ "Retired",
    rowname == "religion2_Other" ~ "Religion (Other)",
    rowname == "religion2_Muslim" ~ "Religion (Muslim)",
    rowname == "religion2_Christian" ~ "Religion (Christian)",
    rowname == "religion2_Atheist/Agnostic" ~ "Religion (Atheist)",
    rowname == "region_QC" ~ "Region (Quebec)",
    rowname == "region_PR" ~ "Region (Prairies)",
    rowname == "region_ON" ~ "Region (Ontario)",
    rowname == "region_BC" ~ "Region (British Columbia)",
    rowname == "region_AT" ~ "Region (Atlantic)",
    rowname == "p_ban_niqab" ~ "Support a ban on the niqab",
    rowname == "income_full_More than $110,000" ~ "Income ($110,000 +)",
    rowname == "income_full_Less than $29,999" ~ "Income (- $29,999)",
    rowname == "income_full_$90,000 - $109,999" ~ "Income ($90,000 - $109,999)",
    rowname == "income_full_$60,000 - $89,999" ~ "Income ($60,000 - $89,999)",
    rowname == "income_full_$30,000 - $59,999" ~ "Income ($30,000 - $59,999)",
    rowname == "french" ~ "French",
    rowname == "female" ~ "Female",
    rowname == "education_No HS" ~ "Education (No High School)",
    rowname == "education_HS" ~ "Education (High School)",
    rowname == "education_Graduate" ~ "Education (Graduate Degree)",
    rowname == "education_Bachelor" ~ "Education (Bachelor Degree)",
    #rowname == "ad_exposure" ~ "Exposure to political advertisements", 
    rowname == "interest" ~ "Interest in the campaign",
    rowname == "knowledge" ~ "Compound knowledge measure",
    TRUE ~ rowname
  )) %>%
  mutate(Diff.Un = abs(Diff.Un)) %>%
  filter(!grepl("<NA>",rowname)) %>%
  ggplot(aes(x = Diff.Un, y = rowname)) +
  geom_point() + 
  theme +
  geom_vline(xintercept = 0.10, linetype = "dashed") +
  labs(x = "Standardised Mean Differences (SMD)", y = "")

ggsave("outputs/figures/figD1.png", width = 7, height = 5)
```

### figD2
```{r}

temp <- filter(CES, date <= "2015-10-19" & !is.na(date), web == 1)

list_of_vars <- c("age","female","french","region","working","retired",
                  "student","religion2","education","income_full",
                  "p_ban_niqab","vote_2011_ndp")

balance_tests <- list()

for (i in 5:10) {
  
  for (j in 5:10) {
    
    if (i <= j) {
      next
    }
    
    #i <- 1
    #j <- 2
    
    group <- paste("group", j, sep = "_")
    
    df <- filter(temp, group %in% c(i, j)) %>%
      .[c(list_of_vars,group)] %>%
      data.frame()
    
    balance_tests[[length(balance_tests)+1]] <- balanceTable(df, treatment = group) %>%
      rownames_to_column(var = "variable")
    
    balance_tests[[length(balance_tests)]]$group <- paste("Week ",i,"-",j, sep = "")
  }

}

balance <- plyr::rbind.fill(balance_tests) %>%
  dplyr::select(group, everything(),-`Fisher/Wilcox Pvalue`) %>%
  rename(p_value = `T-test Pvalue.p.value`,
         smd = `SDiff`,
         t_mean = `Treated Mean`,
         c_mean = `Control Mean`) %>%
  dplyr::filter(!grepl('NA', variable))

mean(abs(balance$smd))
balance %>%
  mutate(balanced = ifelse(abs(smd) <= 0.15,1,0)) %>%
  summarize(prop = sum(balanced)/n())

imbalanced <- balance %>%
  mutate(balanced = ifelse(abs(smd) < 0.10,1,0)) %>%
  group_by(variable, balanced) %>%
  summarize(n = n()) %>%
  spread(key = "balanced", value = "n") %>%
  rename("Imbalanced" = `0`, "Balanced" = `1`) %>%
  mutate(Imbalanced = ifelse(is.na(Imbalanced),0,Imbalanced),
         Balanced = ifelse(is.na(Balanced),0,Balanced)) %>%
  gather(key = "key", value = "value", -variable) %>%
  ungroup() %>%
  mutate(variable = case_when(
    variable == "interest" ~ "Interest in the campaign",
    variable == "knowledge" ~ "Compound knowledge measure",
    variable == "age" ~ "Age",
    variable == "educationBachelor" ~ "Bachelor's Degree",
    variable == "educationGraduate" ~ "Graduate Degree",
    variable == "educationHS" ~ "High School",
    variable == "educationNo HS" ~ "No High School",
    variable == "french" ~ "French",
    variable == "female" ~ "Female",
    variable == "regionAT" ~ "Atlantic Canada",
    variable == "regionBC" ~ "British Columbia",
    variable == "regionON" ~ "Ontario",
    variable == "regionPR" ~ "Prairies",
    variable == "regionQC" ~ "Quebec",
    variable == "income_full$30,000 - $59,999" ~ "$30,000 - $59,999",
    variable == "income_full$60,000 - $89,999" ~ "$60,000 - $89,999",
    variable == "income_full$90,000 - $109,999" ~ "$90,000 - $109,999",
    variable == "income_fullLess than $29,999" ~ "Less than $29,999",
    variable == "income_fullMore than $110,000" ~ "More than $110,000",
    variable == "religion2Atheist/Agnostic" ~ "Atheist/Agnostic",
    variable == "religion2Christian" ~ "Christian",
    variable == "religion2Muslim" ~ "Muslim",
    variable == "religion2Other" ~ "Other Religion",
    variable == "retired" ~ "Retired",
    variable == "student" ~ "Student",
    variable == "working" ~ "Working",
    variable == "p_ban_niqab" ~ "Supports a ban on the niqab",
    variable == "vote_2011_ndp" ~ "Voted for NDP in 2011"
  ))

ggplot(imbalanced, aes(x = variable, y = value, fill = key)) +
  geom_bar(stat = 'identity') +
  theme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y = "Number of week pairs") +
  scale_fill_manual(name = "",values = c("grey","red")) +
  scale_y_continuous(breaks = seq(0,16,2))

ggsave("outputs/figures/figD2.png", width = 7, height = 5)
```

### figD3
```{r}

# Running this will produce numerous The conditional covariance matrix has zero diagonal elements errors. This is because we are doing day-pair comparisons across many variables. Some days do not have people in each cateogry which produces this error (e.g. 52 and 45 Atheist/Agnostic, No HS, $30,000 - $59,999, against the niqab ban in Quebec).

temp <- dplyr::filter(CES, date <= "2015-10-19" & !is.na(date), web == 1)

temp <- dummy_cols(temp, select_columns = "day_of_campaign")

list_of_vars <- c("age","female","french","region","working","retired",
                  "student","religion2","education","income_full",
                  "p_ban_niqab","vote_2011_ndp")

balance_tests <- list()

for (i in 37:69) {
  
  for (j in 37:69) {
    
    if (i <= j) {
      next
    }

    day <- paste("day_of_campaign", j, sep = "_")
    
    df <- filter(temp, day_of_campaign %in% c(i, j)) %>%
      .[c(list_of_vars,day)] %>%
      data.frame()
    
    balance_tests[[length(balance_tests)+1]] <- balanceTable(df, treatment = day) %>%
      rownames_to_column(var = "variable")
    
    balance_tests[[length(balance_tests)]]$day <- paste("Day ",i,"-",j, sep = "")
  }

}

balance <- plyr::rbind.fill(balance_tests) %>%
  dplyr::select(day, everything(),-`Fisher/Wilcox Pvalue`) %>%
  rename(p_value = `T-test Pvalue.p.value`,
         smd = `SDiff`,
         t_mean = `Treated Mean`,
         c_mean = `Control Mean`) %>%
  dplyr::filter(!grepl('NA', variable))

web_daily_mean <- mean(abs(balance$smd))

balance %>%
  mutate(balanced = ifelse(abs(smd) <= 0.15,1,0)) %>%
  summarize(prop = sum(balanced)/n())

imbalanced <- balance %>%
  mutate(balanced = ifelse(abs(smd) <= 0.15,1,0)) %>%
  group_by(variable, balanced) %>%
  summarize(n = n()) %>%
  spread(key = "balanced", value = "n") %>%
  rename("Imbalanced" = `0`, "Balanced" = `1`) %>%
  mutate(Imbalanced = ifelse(is.na(Imbalanced),0,Imbalanced),
         Balanced = ifelse(is.na(Balanced),0,Balanced)) %>%
  gather(key = "key", value = "value", -variable) %>%
  ungroup() %>%
  mutate(variable = case_when(
    variable == "age" ~ "Age",
    variable == "educationBachelor" ~ "Bachelor's Degree",
    variable == "educationGraduate" ~ "Graduate Degree",
    variable == "educationHS" ~ "High School",
    variable == "educationNo HS" ~ "No High School",
    variable == "french" ~ "French",
    variable == "female" ~ "Female",
    variable == "regionAT" ~ "Atlantic Canada",
    variable == "regionBC" ~ "British Columbia",
    variable == "regionON" ~ "Ontario",
    variable == "regionPR" ~ "Prairies",
    variable == "regionQC" ~ "Quebec",
    variable == "income_full$30,000 - $59,999" ~ "$30,000 - $59,999",
    variable == "income_full$60,000 - $89,999" ~ "$60,000 - $89,999",
    variable == "income_full$90,000 - $109,999" ~ "$90,000 - $109,999",
    variable == "income_fullLess than $29,999" ~ "Less than $29,999",
    variable == "income_fullMore than $110,000" ~ "More than $110,000",
    variable == "religion2Atheist/Agnostic" ~ "Atheist/Agnostic",
    variable == "religion2Christian" ~ "Christian",
    variable == "religion2Muslim" ~ "Muslim",
    variable == "religion2Other" ~ "Other Religion",
    variable == "retired" ~ "Retired",
    variable == "student" ~ "Student",
    variable == "working" ~ "Working",
    variable == "p_ban_niqab" ~ "Supports a ban on the niqab",
    variable == "vote_2011_ndp" ~ "Voted for NDP in 2011"
  ))

ggplot(imbalanced, aes(x = variable, y = value, fill = key)) +
  geom_bar(stat = 'identity') +
  theme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y = "Number of day pairs") +
  scale_fill_manual(name = "",values = c("grey","red")) +
  scale_y_continuous(breaks = seq(0,500,50))

ggsave("outputs/figures/figD3.png", width = 7, height = 5)
```

### figD4
```{r}

temp <- filter(CES, date <= "2015-10-19" & !is.na(date), web == 0)

list_of_vars <- c("age","female","french","region","working","retired",
                  "student","religion2","education","income_full",
                  "p_ban_niqab","vote_2011_ndp")

balance_tests <- list()

for (i in 5:10) {
  
  for (j in 5:10) {
    
    if (i <= j) {
      next
    }

    group <- paste("group", j, sep = "_")
    
    df <- filter(temp, group %in% c(i, j)) %>%
      .[c(list_of_vars,group)] %>%
      data.frame()
    
    balance_tests[[length(balance_tests)+1]] <- matchMulti::balanceTable(df, treatment = group) %>%
      rownames_to_column(var = "variable")
    
    balance_tests[[length(balance_tests)]]$group <- paste("Week ",i,"-",j, sep = "")
  }

}

balance <- plyr::rbind.fill(balance_tests) %>%
  dplyr::select(group, everything(),-`Fisher/Wilcox Pvalue`) %>%
  rename(p_value = `T-test Pvalue.p.value`,
         smd = `SDiff`,
         t_mean = `Treated Mean`,
         c_mean = `Control Mean`) %>%
  dplyr::filter(!grepl('NA', variable))

balance %>%
  mutate(balanced = ifelse(abs(smd) <= 0.15,1,0)) %>%
  summarize(prop = sum(balanced)/n())

imbalanced <- balance %>%
  mutate(balanced = ifelse(abs(smd) < 0.10,1,0)) %>%
  group_by(variable, balanced) %>%
  summarize(n = n()) %>%
  spread(key = "balanced", value = "n") %>%
  rename("Imbalanced" = `0`, "Balanced" = `1`) %>%
  mutate(Imbalanced = ifelse(is.na(Imbalanced),0,Imbalanced),
         Balanced = ifelse(is.na(Balanced),0,Balanced)) %>%
  gather(key = "key", value = "value", -variable) %>%
  ungroup() %>%
  mutate(variable = case_when(
    variable == "age" ~ "Age",
    variable == "educationBachelor" ~ "Bachelor's Degree",
    variable == "educationGraduate" ~ "Graduate Degree",
    variable == "educationHS" ~ "High School",
    variable == "educationNo HS" ~ "No High School",
    variable == "french" ~ "French",
    variable == "female" ~ "Female",
    variable == "regionAT" ~ "Atlantic Canada",
    variable == "regionBC" ~ "British Columbia",
    variable == "regionON" ~ "Ontario",
    variable == "regionPR" ~ "Prairies",
    variable == "regionQC" ~ "Quebec",
    variable == "income_full$30,000 - $59,999" ~ "$30,000 - $59,999",
    variable == "income_full$60,000 - $89,999" ~ "$60,000 - $89,999",
    variable == "income_full$90,000 - $109,999" ~ "$90,000 - $109,999",
    variable == "income_fullLess than $29,999" ~ "Less than $29,999",
    variable == "income_fullMore than $110,000" ~ "More than $110,000",
    variable == "religion2Atheist/Agnostic" ~ "Atheist/Agnostic",
    variable == "religion2Christian" ~ "Christian",
    variable == "religion2Muslim" ~ "Muslim",
    variable == "religion2Other" ~ "Other Religion",
    variable == "retired" ~ "Retired",
    variable == "student" ~ "Student",
    variable == "working" ~ "Working",
    variable == "p_ban_niqab" ~ "Supports a ban on the niqab",
    variable == "vote_2011_ndp" ~ "Voted for NDP in 2011"
  ))

ggplot(imbalanced, aes(x = variable, y = value, fill = key)) +
  geom_bar(stat = 'identity') +
  theme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y = "Number of week pairs") +
  scale_fill_manual(name = "",values = c("grey","red")) +
  scale_y_continuous(breaks = seq(0,16,2))

ggsave("outputs/figures/figD4.png", width = 7, height = 5)
```

### figD5
```{r}

temp <- dplyr::filter(CES, date <= "2015-10-19" & !is.na(date), web == 0)

temp <- dummy_cols(temp, select_columns = "day_of_campaign")

list_of_vars <- c("age","female","french","region","working","retired",
                  "student","religion2","education","income_full",
                  "p_ban_niqab","vote_2011_ndp")

balance_tests <- list()

for (i in 37:69) {
  
  for (j in 37:69) {
    
    if (i <= j) {
      next
    }
  
    day <- paste("day_of_campaign", j, sep = "_")
    
    df <- filter(temp, day_of_campaign %in% c(i, j)) %>%
      .[c(list_of_vars,day)] %>%
      data.frame()
    
    balance_tests[[length(balance_tests)+1]] <- balanceTable(df, treatment = day) %>%
      rownames_to_column(var = "variable")
    
    balance_tests[[length(balance_tests)]]$day <- paste("Day ",i,"-",j, sep = "")
  }

}

balance <- plyr::rbind.fill(balance_tests) %>%
  dplyr::select(day, everything(),-`Fisher/Wilcox Pvalue`) %>%
  rename(p_value = `T-test Pvalue.p.value`,
         smd = `SDiff`,
         t_mean = `Treated Mean`,
         c_mean = `Control Mean`) %>%
  dplyr::filter(!grepl('NA', variable))

web_daily_mean <- mean(abs(balance$smd))

balance %>%
  mutate(balanced = ifelse(abs(smd) <= 0.15,1,0)) %>%
  summarize(prop = sum(balanced)/n())

imbalanced <- balance %>%
  mutate(balanced = ifelse(abs(smd) <= 0.15,1,0)) %>%
  group_by(variable, balanced) %>%
  summarize(n = n()) %>%
  spread(key = "balanced", value = "n") %>%
  rename("Imbalanced" = `0`, "Balanced" = `1`) %>%
  mutate(Imbalanced = ifelse(is.na(Imbalanced),0,Imbalanced),
         Balanced = ifelse(is.na(Balanced),0,Balanced)) %>%
  gather(key = "key", value = "value", -variable) %>%
  ungroup() %>%
  mutate(variable = case_when(
    variable == "age" ~ "Age",
    variable == "educationBachelor" ~ "Bachelor's Degree",
    variable == "educationGraduate" ~ "Graduate Degree",
    variable == "educationHS" ~ "High School",
    variable == "educationNo HS" ~ "No High School",
    variable == "french" ~ "French",
    variable == "female" ~ "Female",
    variable == "regionAT" ~ "Atlantic Canada",
    variable == "regionBC" ~ "British Columbia",
    variable == "regionON" ~ "Ontario",
    variable == "regionPR" ~ "Prairies",
    variable == "regionQC" ~ "Quebec",
    variable == "income_full$30,000 - $59,999" ~ "$30,000 - $59,999",
    variable == "income_full$60,000 - $89,999" ~ "$60,000 - $89,999",
    variable == "income_full$90,000 - $109,999" ~ "$90,000 - $109,999",
    variable == "income_fullLess than $29,999" ~ "Less than $29,999",
    variable == "income_fullMore than $110,000" ~ "More than $110,000",
    variable == "religion2Atheist/Agnostic" ~ "Atheist/Agnostic",
    variable == "religion2Christian" ~ "Christian",
    variable == "religion2Muslim" ~ "Muslim",
    variable == "religion2Other" ~ "Other Religion",
    variable == "retired" ~ "Retired",
    variable == "student" ~ "Student",
    variable == "working" ~ "Working",
    variable == "p_ban_niqab" ~ "Supports a ban on the niqab",
    variable == "vote_2011_ndp" ~ "Voted for NDP in 2011"
  ))

ggplot(imbalanced, aes(x = variable, y = value, fill = key)) +
  geom_bar(stat = 'identity') +
  theme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y = "Number of day pairs") +
  scale_fill_manual(name = "",values = c("grey","red")) +
  scale_y_continuous(breaks = seq(0,500,50))

ggsave("outputs/figures/figD5.png", width = 7, height = 5)

```

## D-3 Pre-period linear trend and placebo test

### tabD1
```{r}

CES <- mutate(CES, court_linear_alt = as.numeric(date - 16692))

pre_qc <- list()

pre_qc[[1]] <- lm(data = filter(CES, web == 1, date <= "2015-09-15"), 
                100*ndp_vote ~ female + age + french + 
                  region_AT + region_PR + region_BC +working + student + retired + 
                  education+ vote_2011_ndp + court_linear_alt*quebec)
pre_qc[[2]] <- lm(data = filter(CES, web == 1, date <= "2015-09-15"), 
                100*ndp_vote ~ court_linear_alt*quebec)

screenreg(pre_qc, single.row = TRUE)

texreg(pre_qc, single.row = TRUE,
       scalebox = 0.7, float.pos = "ht!", label = "tabD1", dcolumn = TRUE,  use.packages = FALSE,
       caption = "Pre-trend test for DID comparing Quebec and the Rest of Canada", caption.above = TRUE,
       custom.note = ("\\parbox{1\\linewidth}{\\vspace{2pt}%stars. OLS estimations. Standard errors in parentheses. Dependent variable is vote for the NDP in the 2015 Canadian federal election (binary).}"),
       custom.model.names = c("1: Pre-trend test (controls)",
                              "2: Pre-trend test (no controls)"),
       custom.coef.names = c("Constant","Female","Age","French",
                             "British Columbia","Ontario","Prairies","Quebec",
                             "Working","Student","Retired",
                             "High School","Bachelor's Degree","Graduate Degree",
                             "Vote 2011 NDP","Pre-trend","Pre-trend:Quebec"),
       file = "outputs/tables/tabD1.tex")

```

### tabD2
```{r}

CES <- mutate(CES,
               court_linear_alt = as.numeric(date - 16692),
               pre_trend = ifelse(date > "2015-09-08", 1, 0))

pre_qc_1 <- list()

pre_qc_1[[1]] <- lm(data = filter(CES, web == 1, date <= "2015-09-15"), 
                100*ndp_vote ~ court_linear_alt*quebec)
pre_qc_1[[2]] <- lm(data = filter(CES, web == 1, date <= "2015-09-15"), 
                100*ndp_vote ~ pre_trend*quebec)
pre_qc_1[[3]] <- lm(data = filter(CES, web == 1, date <= "2015-09-15"), 
                100*ndp_vote ~ female + age + french + region_AT + region_PR + region_BC + quebec + working + student + retired + 
                  education+ vote_2011_ndp + pre_trend*quebec)

screenreg(pre_qc_1, single.row = TRUE)

texreg(pre_qc_1, single.row = TRUE, dcolumn = TRUE,  use.packages = FALSE,
       scalebox = 0.7, float.pos = "ht!", label = "tabD2",
       caption = "Pre-trend test for DID comparing Quebec and the Rest of Canada", caption.above = TRUE,
       custom.note = ("\\parbox{1\\linewidth}{\\vspace{2pt}%stars. OLS estimations. Standard errors in parentheses. Dependent variable is vote for the NDP in the 2015 Canadian federal election (binary). Sub-sample of respondents before September 15.}"),
       custom.model.names = c("1: Trend",
                              "2: 1-week dummy",
                              "3: 1-week dummy (controls)"),
       custom.coef.names = c("Constant",
                             "Pre-trend","Quebec","Pre-trend x Quebec",
                             "Sep 8","Sep 8 x Quebec",
                             "Female","Age","French",
                             "Atlantic","Prairies",
                             "British Columbia",
                             "Working","Student","Retired",
                             "High School","Bachelor's Degree","Graduate Degree",
                             "Vote 2011 NDP","Pre-trend x Quebec"),
       file = "outputs/tables/tabD2.tex")

```

## E-1 Regression discontinuity

### tabE1
```{r}

CES <- mutate(CES, court_linear_alt = as.numeric(date - 16693))

for (k in c("Quebec")) {
  
  df <- data.frame(sample = NA,
                 bw_selection = NA,
                 bw = NA,
                 p = NA,
                 coef = NA,
                 se = NA,
                 pv = NA)

  if (k == "Quebec") {
    
    x <- filter(CES, quebec == 1) %>% pull(court_linear_alt)
    y <- filter(CES, quebec == 1) %>% pull(ndp_vote)
    name = "qc"

  } 
  
  for (j in c(1:3)) { 
    
    for (i in c(3,5,7,9)) {
      
      
      if (i == 3 & j != 1 | i == 5 & j == 3) {
        next
      }
      
      rd_model <- rdrobust(y = y, x = x, 
                           kernel = "triangular", scaleregul = 1,
                           p = j,
                           h = i)
      
      df[nrow(df)+1,"coef"] <- rd_model$coef[3]
      df[nrow(df),"bw_selection"] <- "manual"
      df[nrow(df),"se"] <- rd_model$se[3]
      df[nrow(df),"pv"] <- rd_model$pv[3]
      df[nrow(df),"p"] <- j
      df[nrow(df),"bw"] <- i
      df[nrow(df),"sample"] <- k

    }
    
    rd_model <- rdrobust(y = y, x = x, 
                         kernel = "triangular", scaleregul = 1,
                         p = j,
                         bwselect = "mserd")
    
    df[nrow(df)+1,"coef"] <- rd_model$coef[3]
    df[nrow(df),"bw_selection"] <- "mserd"
    df[nrow(df),"se"] <- rd_model$se[3]
    df[nrow(df),"pv"] <- rd_model$pv[3]
    df[nrow(df),"p"] <- j
    df[nrow(df),"bw"] <- round(rd_model$bws[1,1],2)
    df[nrow(df),"sample"] <- k
    
  }
  
  df <- df[-1,]
  
  colnames(df) <- c("Sub-sample","Selection","Bandwidth","Poly. order","Coefficient","St. Err.","P-value")
  
  print(xtable::xtable(df, caption = "Testing RD around September 15th with a variety polynomials of bin widths", label = "tabE1"),
        file= "outputs/tables/tabE1.tex",
        scalebox=0.9, caption.placement = "top",
        include.rownames=FALSE)
  
}
```

### figE1/2
```{r}

x <- filter(CES, quebec == 1) %>% pull(court_linear_alt)
y <- filter(CES, quebec == 1) %>% pull(ndp_vote)
h <- 12.58
name <- "qc"

png(filename = "outputs/figures/figE1.png",
    res = 72, width = 6, height = 4, units = "in")
rdplot(y = y, x = x, p = 1, kernel = "triangular", x.lim = c(-10, 10), h = 7,
       x.label = "Days surrounding September 15", y.label = "Mean vote intention for NDP", title = "")
dev.off()

png(filename = "outputs/figures/figE2.png",
    res = 72, width = 6, height = 4, units = "in")
rdplot(y = y, x = x, p = 1, kernel = "triangular", x.lim = c(-h-1,h+1),
       x.label = "Days surrounding September 15", y.label = "Mean vote intention for NDP", title = "")
dev.off()


```

## E-2 DID 7-day window


### tabE2
```{r}

models <- list()
model_names <- vector()

for (i in as.Date("2015-09-01"):as.Date("2015-10-08")) {

  # Only do every third day
  j <- as.numeric(as.Date(i) - as.Date("2015-08-31"))
  if (j %% 3 != 0) next 
    w <- j/3
  
  model_names[w] <- i
  
  # Create a temporary dataframe with 14 day window on either side of
  # a fictitious treatment and add a categorical and linear trend from
  # that date
  df_temp <- filter(CES, web == 1,
                    date > (as.Date(i)-7),
                    date <= (as.Date(i)+7)) %>%
    mutate(treatment = ifelse(date > as.Date(i), 1, 0),
           treatment_linear = ifelse(treatment != 1, 0,
                                     as.numeric(date - as.Date(i))))

  # Run a model on the data subset
  models[[w]] <- lm(data = df_temp, 
                    100*ndp_vote ~ treatment*quebec)
    
}

model_names <- as.character(as.Date(model_names))

texreg(models, custom.model.names = model_names, dcolumn = TRUE,  use.packages = FALSE,
       scalebox = 0.6, float.pos = "ht!", caption.above = TRUE, label = "tabE2",
       caption = "Testing alternative treatment dates with pre and post windows of 7 days (Quebec versus Rest of Canada) and no controls",
              custom.note = ("\\parbox{1.3\\linewidth}{\\vspace{2pt}%stars. OLS estimations (DID models). Standard errors in parentheses. Dependent variable is vote for the NDP in the 2015 Canadian federal election (binary). The name of the model shows the date around which the 7-day pre- and post-window is centered.}"),
        custom.coef.names = c("Constant","Treatment","Quebec","Treatment:Quebec"),
        file = "outputs/tables/tabE2.tex")
```



### tabE3
```{r}


models <- list()
model_names <- vector()

for (i in as.Date("2015-09-01"):as.Date("2015-10-08")) {

  # Only do every third day
  j <- as.numeric(as.Date(i) - as.Date("2015-08-31"))
  if (j %% 3 != 0) next 
    w <- j/3
  
  model_names[w] <- i
  
  # Create a temporary dataframe with 14 day window on either side of
  # a fictitious treatment and add a categorical and linear trend from
  # that date
  df_temp <- filter(CES, web == 1,
                    date > (as.Date(i)-7),
                    date <= (as.Date(i)+7)) %>%
    mutate(treatment = ifelse(date > as.Date(i), 1, 0),
           treatment_linear = ifelse(treatment != 1, 0,
                                     as.numeric(date - as.Date(i))))

  # Run a model on the data subset
  models[[w]] <- lm(data = df_temp, 
                    100*ndp_vote ~ female + age + french + 
                      region_AT + region_PR + region_BC + working + student + retired + education + 
                      vote_2011_ndp + treatment*quebec)
    
}

model_names <- as.character(as.Date(model_names))

texreg(models, custom.model.names = model_names, dcolumn = TRUE,  use.packages = FALSE,
       scalebox = 0.6, float.pos = "ht!", caption.above = TRUE, label = "tabE3",
       caption = "Testing alternative treatment dates with pre and post windows of 7 days (Quebec versus Rest of Canada)",
              custom.note = ("\\parbox{1.3\\linewidth}{\\vspace{2pt}%stars. OLS estimations (DID models). Standard errors in parentheses. Dependent variable is vote for the NDP in the 2015 Canadian federal election (binary). The name of the model shows the date around which the 7-day pre- and post-window is centered.}"),
        custom.coef.names = c("Constant","Female","Age",
                              "French","Atlantic","Prairies","British Columbia",
                              "Full-time worker","Student","Retired",
                              "High School","Bachelor","Graduate",
                              "Vote NDP 2011", "Treatment", "Quebec", "Treatment:Quebec"),
        file = "outputs/tables/tabE3.tex")

```



## F-1 Sensitivity analysis

### konfound
```{r}

filter(CES, web == 1) %>%
  mutate(court_quebec = court*quebec) %>%
  lm(data = ., ndp_vote ~ female + age + french + 
                  region_AT + region_PR + region_BC +working + student + retired + education +
                  vote_2011_ndp + court_quebec) %>%
  konfound(., court_quebec)

filter(CES, web == 1) %>%
  mutate(court_linear_quebec = court_linear*quebec) %>%
  lm(data = ., ndp_vote ~ female + age + french + 
                  region_AT + region_PR + region_BC +working + student + retired + education +
                  vote_2011_ndp + court_linear_quebec) %>%
  konfound(., court_linear_quebec)

```

## F-2 Logistic Regression

### tabF1
```{r}

tab <- list()

tab[[1]] <- glm(data = filter(CES, web == 1), family = "binomial",
                       ndp_vote ~ female + french + age + region_AT + region_PR + region_BC +
                         working + student + retired + education +
                         vote_2011_ndp + court*quebec)

tab[[2]] <- glm(data = filter(CES, web == 1),  family = "binomial",
                       ndp_vote ~ female + age + french + 
                         region_AT + region_PR + region_BC +working + student + retired + education + 
                         vote_2011_ndp + court*quebec + court_linear*quebec)

tab[[3]] <-  glm.cluster(data = filter(CES, web == 1), 
                         cluster = "date",family="binomial",
                         formula = ndp_vote ~ female + age + french + 
                           region_AT + region_PR + region_BC +working + student + retired + education +
                           vote_2011_ndp + niqab_7*quebec)

screenreg(list(tab[[1]],tab[[2]],tab[[3]]$glm_res), include.ci = FALSE,
          omit.coef = "(female)|(age)|(french)|(region)|
          (PR)|(working)|(student)|(retired)|(education)")

texreg(list(tab[[1]],tab[[2]],tab[[3]]$glm_res),
       scalebox = 1, float.pos = "ht!", label = "tabF1", include.ci = FALSE,
       dcolumn = TRUE,  use.packages = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP", 
       caption.above = TRUE, single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)",
       custom.note = ("\\parbox{1.0\\linewidth}{\\vspace{2pt}%stars. Linear probability models for DID estimations with robust standard errors for Models 1 and 2 and clustered standard errors for Model 3 in parentheses. Dependent variable is vote intention for the NDP (binary variable). All models use full web sample.}"), 
       custom.model.names = c("1: Binary DID","2: Linear Trend ","3: 7-day media"),
       custom.coef.names = c("Constant","Quebec","Voted NDP 2011",
                             "Ruling","Ruling x Quebec","Trend","Trend x Quebec",
                             "7-day niqab","7-day niqab x Quebec"),
       reorder.coef = c(1,2,3,4,6,8,5,7,9),
       groups = list(" " = 1:6,"DID effects" = 7:9),
       file = "outputs/tables/tabF1.tex")

```

## F-3 Alternative samples

### tabF2
```{r}

tab1 <- list()

tab1[[1]] <- lm_robust(data = filter(CES) %>% mutate(ndp_vote = ndp_vote*100),
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + french + age + region_AT + region_PR + region_BC +web + 
                         working + student + retired + education +
                         vote_2011_ndp + court*quebec)

tab1[[2]] <- lm_robust(data = filter(CES) %>% mutate(ndp_vote = ndp_vote*100), 
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + age + french + web + 
                         region_AT + region_PR + region_BC +working + student + retired + education + 
                         vote_2011_ndp + court*quebec + court_linear*quebec)

tab1[[3]] <-  lm.cluster(data = filter(CES) %>% mutate(ndp_vote = ndp_vote*100), 
                         cluster = "date",
                         formula = ndp_vote ~ female + age + french + web + 
                           region_AT + region_PR + region_BC +working + student + retired + education +
                           vote_2011_ndp + niqab_7*quebec)

screenreg(list(tab1[[1]],tab1[[2]],tab1[[3]]$lm_res), include.ci = FALSE,
          omit.coef = "(female)|(age)|(french)|(region)|
          (PR)|(working)|(student)|(retired)|(education)")

texreg(list(tab1[[1]],tab1[[2]],tab1[[3]]$lm_res),
       scalebox = 1, float.pos = "ht!", label = "tabF2", include.ci = FALSE,
       dcolumn = TRUE,  use.packages = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP", 
       caption.above = TRUE, single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)",
       custom.note = ("\\parbox{1.0\\linewidth}{\\vspace{2pt}%stars. Linear probability models for DID estimations with robust standard errors for Models 1 and 2 and clustered standard errors for Model 3 in parentheses. Dependent variable is vote intention for the NDP (binary variable). All models use full web sample.}"), 
       custom.model.names = c("1: Binary DID","2: Linear Trend ","3: 7-day media"),
       custom.coef.names = c("Constant","Web","Voted NDP 2011",
                             "Ruling","Quebec","Ruling x Quebec","Trend","Trend x Quebec",
                             "7-day niqab","7-day niqab x Quebec"),
       reorder.coef = c(1,2,3,4,5,7,9,6,8,10),
       groups = list(" " = 1:7,"DID effects" = 8:10),
       file = "outputs/tables/tabF2.tex")
```

## F-4 Local Parliament Project Data

### tabF3
```{r}

tab_lpp <- list()

tab_lpp[[1]] <- lm_robust(data = LPP %>% mutate(ndp_vote = ndp_vote*100),
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + french + age + region_AT + region_PR + region_BC +
                         working + student + retired + education +
                         vote_2011_ndp + court*quebec)

tab_lpp[[2]] <- lm_robust(data = LPP %>% mutate(ndp_vote = ndp_vote*100), 
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + age + french + 
                         region_AT + region_PR + region_BC +working + student + retired + education + 
                         vote_2011_ndp + court*quebec + court_linear*quebec)

tab_lpp[[3]] <-  lm.cluster(data = LPP %>% mutate(ndp_vote = ndp_vote*100), 
                         cluster = "date",
                         formula = ndp_vote ~ female + age + french + 
                           region_AT + region_PR + region_BC +working + student + retired + education +
                           vote_2011_ndp + niqab_7*quebec)

screenreg(list(tab_lpp[[1]],tab_lpp[[2]],tab_lpp[[3]]$lm_res), include.ci = FALSE,
          omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)")

texreg(list(tab_lpp[[1]],tab_lpp[[2]],tab_lpp[[3]]$lm_res), include.ci = FALSE,
       scalebox = 1, float.pos = "ht!", label = "tabF3",
       dcolumn = TRUE,  use.packages = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP", 
       caption.above = TRUE, single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)",
       custom.note = ("\\parbox{1.0\\linewidth}{\\vspace{2pt}%stars. Linear probability models for DID estimations with robust standard errors for Models 1 and 2 and clustered standard errors for Model 3 in parentheses. Dependent variable is vote intention for the NDP (binary variable). All models use full LPP sample.}"), 
       custom.model.names = c("1: Binary DID","2: Linear Trend ","3: 7-day media"),
       custom.coef.names = c("Constant","Voted NDP 2011",
                             "Ruling","Quebec","Ruling x Quebec","Trend","Trend x Quebec",
                             "7-day niqab","7-day niqab x Quebec"),
       reorder.coef = c(1,2,3,4,6,8,5,7,9),
       groups = list(" " = 1:6,"DID effects" = 7:9),
       file = "outputs/tables/tabF3.tex")

```

## F-5 Quebec effect heterogeneity

### tabF4
```{r}

CES_exposure = filter(CES, web == 1, quebec == 1) %>% mutate(ndp_vote = ndp_vote*100)

tabs = list()

tabs[[1]] = lm_robust(data = CES_exposure,
               ndp_vote ~ female + age + french + 
                 working + student + retired + education +
                 vote_2011_ndp + court*p_media_exposure)
tabs[[2]] = lm_robust(data = CES_exposure,
               ndp_vote ~ female + age + french + 
                 working + student + retired + education +
                 vote_2011_ndp + court*knowledge)
tabs[[3]] = lm_robust(data = CES_exposure,
               ndp_vote ~ female + age + french + 
                 working + student + retired + education +
                 vote_2011_ndp + isslist_immg)

screenreg(tabs, include.ci = FALSE,
          omit.coef = "(female)|(age)|(french)|(ON)|(BC)|(PR)|(PR)|(working)|(student)|(retired)|(education)")

texreg(tabs, include.ci = FALSE,
       scalebox = 0.80, float.pos = "ht!", label = "tabF4",
       dcolumn = TRUE,  use.packages = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP, moderated by various measures of exposure to the niqab story", 
       caption.above = TRUE, #single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(ON)|(BC)|(PR)|(PR)|(working)|(student)|(retired)|(education)",
       custom.note = ("\\parbox{1.25\\linewidth}{\\vspace{2pt}%stars. Linear probability models for DID estimations with standard errors in parentheses. Dependent variable is vote intention for the NDP in the 2015 Canadian federal election (binary). All models use full Quebec sample.}"), 
       custom.model.names = c("1: Media Exposure",
                              "2: Political knowledge",
                              "3: Immigration-news"),
       custom.coef.names = c("Constant",
                             "Voted NDP 2011",
                             "Ruling",
                             "Media exposure",
                             "Ruling: Media exposure",
                             "Knowledge",
                             "Ruling x Knowledge",
                             "Heard about immigration in past week"),
       file = "outputs/tables/tabF4.tex")
```

### tabF5
```{r}

exp_descriptives <- filter(CES, web == 1, quebec == 1) %>% 
  dplyr::select(p_media_exposure,knowledge, isslist_immg) %>%
  psych::describe(skew = TRUE, check = TRUE) %>%
  dplyr::select(n, mean, sd, median, min, max) %>%
  as_tibble() %>%
  rownames_to_column(.)

colnames(exp_descriptives) <- c("Variable","n","Mean","SD","Median","Min","Max")

exp_descriptives[,'Variable'] <- c("1: Post-Election Media Exposure",
                                   "2: Knowledge battery",
                                   "3: Exposure to immigration related news")

exp_descriptives <- arrange(exp_descriptives, Variable)

print(xtable::xtable(exp_descriptives, 
                     digits = c(0,0,0,2,2,2,0,0),label = "tabF5",
                     caption = "Descriptive statistics for exposure measures"),
      file="outputs/tables/tabF5.tex", 
      caption.placement = "top", scalebox='0.9', include.rownames = FALSE)
```

## G-1 Visualizations of media tone and mentions

### figG1
```{r}

daily_mentions = media_sentiment %>% 
  group_by(date) %>% 
  summarize(ndp_mentions = mean(ndp_mentions), lpc_mentions = mean(lpc_mentions))

ggplot(daily_mentions, aes(x = date)) +
  geom_ma(aes(y = ndp_mentions, col = "NDP-associated terms"), ma_fun = SMA, n = 7, lty = 1, cex = 1) +
  geom_ma(aes(y = lpc_mentions, col = "Liberal-associated terms"), ma_fun = SMA, n = 7, lty = 1, cex = 1) +
  scale_x_date(breaks = seq(as.Date("2015-09-15"),as.Date("2015-10-19"),7),
               limits = c(as.Date("2015-09-08"),as.Date("2015-10-19")),
               date_labels = "%b %e") +
  theme + 
  scale_color_manual(name = "",values = c("#D71920","#F37021")) +
  labs(x = "", y = "Average number of mentions per story") + 
  geom_hline(aes(yintercept = mean(daily_mentions$ndp_mentions)),
             linetype = 2, size = 1) + 
  annotate("text", x = as.Date("2015-09-12"),
           y = mean(daily_mentions$ndp_mentions)+0.5,
           col = "black",
           label = "Daily mean NDP mentions (3.2)", size = 2) + 
  geom_hline(aes(yintercept = mean(daily_mentions$lpc_mentions)),
             linetype = 2, size = 1) +
  annotate("text", x = as.Date("2015-09-12"),
           y = mean(daily_mentions$lpc_mentions)+0.5,
           col = "black",
           label = "Daily mean Liberal mentions (1.9)", size = 2)

ggsave("outputs/figures/figG1.png", width = 5, height = 3)
```

### figG2
```{r}

ggplot(daily_sentiment, aes(x = date)) +
  geom_ma(aes(y = prop_neg_lpc, col = "Liberal-associated mentions"),
          size = 1, linetype = "solid", n = 7) +
  geom_ma(aes(y = prop_neg_ndp, col = "NDP-associated mentions"),
          size = 1, linetype = "solid", n = 7) +
  scale_x_date(breaks = seq(as.Date("2015-09-15"), as.Date("2015-10-19"),7),
               limits = c(as.Date("2015-09-15"),as.Date("2015-10-19")),
               date_labels = "%b %e") +
  theme + 
  scale_color_manual(name = "",values = c("#D71920","#F37021")) +
  labs(x = "",
       y = "Proportion of negative sentiment \n directed towards Party or Party Leaders") +
  geom_hline(aes(yintercept = mean(daily_sentiment$prop_neg_ndp, na.rm = TRUE)),
             linetype = 2, size = 0.5) + 
  annotate("text", x = as.Date("2015-09-19"),
           y = mean(daily_sentiment$prop_neg_ndp, na.rm = TRUE)+0.01,
           col = "black",
           label = "Mean NDP sentiment (0.51)", size = 2) + 
  geom_hline(aes(yintercept = mean(daily_sentiment$prop_neg_lpc, na.rm = TRUE)),
             linetype = 2, size = 0.5) +
  annotate("text", x = as.Date("2015-09-19"),
           y = mean(daily_sentiment$prop_neg_lpc, na.rm = TRUE)+0.01,
           col = "black",
           label = "Mean Liberal sentiment (0.46)", size = 2)

ggsave("outputs/figures/figG2.png", width = 5, height = 3)
```

## H-1 The natural experiment modelling strategy

### tabH1
```{r}

alt <- list()

alt[[1]] <- lm(data = filter(CES, web == 1), 
   100*ndp_vote ~ court)

alt[[2]] <- lm(data = filter(CES, web == 1), 
   100*ndp_vote ~ court + female + age + french + 
                  region_AT + region_PR + region_BC +working + student + retired + education +
                  vote_2011_ndp)

alt[[3]] <- lm(data = filter(CES, web == 1), 
   100*ndp_vote ~ court*quebec)

alt[[4]] <- lm(data = filter(CES, web == 1), 
   100*ndp_vote ~ court*quebec + female + age + french + 
                  region_AT + region_PR + region_BC + working + student + retired + education +
                  vote_2011_ndp)

texreg(alt, scalebox = 0.8, float.pos = "ht!", label = "tabH1", 
       dcolumn = TRUE,  use.packages = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP", 
       caption.above = TRUE, single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(region)|(working)|(student)|(retired)|(education)|(vote_2011_ndp)",
       custom.note = ("\\parbox{1.0\\linewidth}{\\vspace{2pt}%stars. Linear probability models based on a natural experiment setup with dependent variable as vote intention for the NDP in the 2015 Canadian federal election (binary). All models use full web sample.}"),
       custom.model.names = c("1: Overall","2: Overall",
                              "3: Quebec", "4: Quebec"),
       custom.coef.names = c("Constant",
                             "Ruling","Quebec",
                             "Ruling x Quebec"),
       reorder.coef = c(1,3,2,4),
       file = "outputs/tables/tabH1.tex")
```

## H-2 Natural experiment modelling strategy - panel data

### descriptive
```{r}
CES %>%
  filter(!is.na(p_ban_niqab)) %>% 
  group_by(quebec, p_ban_niqab, court) %>%
  summarize(ndp_vote = mean(ndp_vote, na.rm = TRUE),
            ndp_voted = mean(ndp_voted, na.rm = TRUE),
            b = n())
```

### figH1
```{r}
rbind(
  t.test(filter(CES, quebec == 1, court == 0, p_ban_niqab == 0)$ndp_voted,
         filter(CES, quebec == 1, court == 0, p_ban_niqab == 0)$ndp_vote) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 0, p_ban_niqab = 0),
  t.test(filter(CES, quebec == 1, court == 1, p_ban_niqab == 0)$ndp_voted,
         filter(CES, quebec == 1, court == 1, p_ban_niqab == 0)$ndp_vote) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 1, p_ban_niqab = 0),
  t.test(filter(CES, quebec == 1, court == 0, p_ban_niqab == 1)$ndp_voted,
         filter(CES, quebec == 1, court == 0, p_ban_niqab == 1)$ndp_vote) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 0, p_ban_niqab = 1),
  t.test(filter(CES, quebec == 1, court == 1, p_ban_niqab == 1)$ndp_voted,
         filter(CES, quebec == 1, court == 1, p_ban_niqab == 1)$ndp_vote) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 1, p_ban_niqab = 1)
) %>%
  mutate(p_ban_niqab = ordered(ifelse(p_ban_niqab == 1, "and in FAVOUR of the ban","and AGAINST the ban"),
                               levels = c("and in FAVOUR of the ban","and AGAINST the ban")),
         court = factor(ifelse(court == 1,
                                 "Took survey AFTER",
                                 "Took survey BEFORE"),
                          ordered = TRUE, 
                          levels = c("Took survey BEFORE",
                                     "Took survey AFTER")),
         text = paste(court, p_ban_niqab, sep = "\n")) %>%
  ggplot(aes(x = text, y = estimate1-estimate2)) +
  geom_point() + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) + 
  theme_bw() +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(x = "", y = "Percentage point change in vote support\nfor NDP in Quebec")

ggsave("outputs/figures/figH1.png", width = 5, height = 2.5)

```

### figH2
```{r}

rbind(
  t.test(filter(CES, quebec == 1, court == 0, p_ban_niqab == 0)$p_ndp_pid,
         filter(CES, quebec == 1, court == 0, p_ban_niqab == 0)$ndp_pid) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 0, p_ban_niqab = 0),
  t.test(filter(CES, quebec == 1, court == 1, p_ban_niqab == 0)$p_ndp_pid,
         filter(CES, quebec == 1, court == 1, p_ban_niqab == 0)$ndp_pid) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 1, p_ban_niqab = 0),
  t.test(filter(CES, quebec == 1, court == 0, p_ban_niqab == 1)$p_ndp_pid,
         filter(CES, quebec == 1, court == 0, p_ban_niqab == 1)$ndp_pid) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 0, p_ban_niqab = 1),
  t.test(filter(CES, quebec == 1, court == 1, p_ban_niqab == 1)$p_ndp_pid,
         filter(CES, quebec == 1, court == 1, p_ban_niqab == 1)$ndp_pid) %>%
    tidy() %>% dplyr::select(estimate1,estimate2,conf.low,conf.high) %>%
    mutate(court = 1, p_ban_niqab = 1)
) %>%
  mutate(p_ban_niqab = ordered(ifelse(p_ban_niqab == 1, "and in FAVOUR of the ban","and AGAINST the ban"),
                               levels = c("and in FAVOUR of the ban","and AGAINST the ban")),
         court = factor(ifelse(court == 1,
                                 "Took survey AFTER",
                                 "Took survey BEFORE"),
                          ordered = TRUE, 
                          levels = c("Took survey BEFORE",
                                     "Took survey AFTER")),
         text = paste(court, p_ban_niqab, sep = "\n")) %>%
  ggplot(aes(x = text, y = estimate1-estimate2)) +
  geom_point() + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) + 
  theme_bw() +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(x = "", y = "Percentage point change in NDP\npartisanship in Quebec")

ggsave("outputs/figures/figH2.png", width = 5, height = 2.5)

```

## H-3 NDP 2011 Voters

### tabH2
```{r}

tab_ndp_dat = filter(CES, web == 1, quebec == 1) %>% mutate(ndp_vote = ndp_vote*100)

tab_ndp <- list()

tab_ndp[[1]] <- lm_robust(data = tab_ndp_dat,
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + french + age +
                         working + student + retired + education +
                         vote_2011_ndp + court*vote_2011_ndp)

tab_ndp[[2]] <- lm_robust(data = tab_ndp_dat, 
                       ci = TRUE, se_type = "HC2",
                       ndp_vote ~ female + age + french + 
                         working + student + retired + education + 
                         vote_2011_ndp + court*vote_2011_ndp + court_linear*vote_2011_ndp)

tab_ndp[[3]] <-  lm.cluster(data = tab_ndp_dat, 
                         cluster = "date",
                         formula = ndp_vote ~ female + age + french + 
                           working + student + retired + education +
                           vote_2011_ndp + niqab_7*vote_2011_ndp)

screenreg(list(tab_ndp[[1]],tab_ndp[[2]],tab_ndp[[3]]$lm_res), include.ci = FALSE,
          omit.coef = "(female)|(age)|(french)|(ON)|(BC)|(PR)|
          (PR)|(working)|(student)|(retired)|(education)")

texreg(list(tab_ndp[[1]],tab_ndp[[2]],tab_ndp[[3]]$lm_res),
       scalebox = 0.8, float.pos = "ht!", label = "tabH2", include.ci = FALSE,
       dcolumn = TRUE,  use.packages = FALSE,
       caption = "The effects of the niqab ruling on vote intention for the NDP", 
       caption.above = TRUE, single.row = TRUE,
       omit.coef = "(female)|(age)|(french)|(ON)|(BC)|(PR)|(PR)|(working)|(student)|(retired)|(education)",
       custom.note = ("\\parbox{1.15\\linewidth}{\\vspace{2pt}%stars. Linear probability models for DID estimations with robust standard errors for Models 1 and 2 and clustered standard errors for Model 3 in parentheses. Dependent variable is vote intention for the NDP (binary variable). All models use full web sample.}"), 
       custom.model.names = c("1: Binary DID","2: Linear Trend ","3: 7-day media"),
       custom.coef.names = c("Constant","Voted NDP 2011",
                             "Ruling","Ruling x Voted NDP 2011","Trend","Trend x Voted NDP 2011",
                             "7-day niqab","7-day niqab x Voted NDP 2011"),
       reorder.coef = c(1,2,3,5,7,4,6,8),
       groups = list(" " = 1:5,"DID effects" = 6:8),
       file = "outputs/tables/tabH2.tex")

```
